#include <iostream>
#include <vector>
#include <cmath>
#include <sstream>
#include "../Eigen/Dense"
#include "../spline_functions.h"
// #include "matplotlibcpp.h"

// namespace plt = matplotlibcpp;

using namespace std;

// 定义函数 f(x)
double f(double x) {
    return 1.0 / (1 + x * x);
}

// 计算插值误差的函数
double calculateError(double x, const BSpline& spline) {
    return fabs(spline.evaluate(x) - f(x));
}

int main() {
    std::cout << "-------Problem D: Testing the interpolation error of the two cardinal B-spline interpolants-------" << std::endl;
    // 定义节点
    std::vector<double> x_values_1, x_values_2;
    std::vector<double> y_values_1, y_values_2;

    // 为定理3.57（三次B样条）定义节点
    int N1 = 11;
    for (int i = 1; i <= N1; ++i) {
        x_values_1.push_back(-6 + i);
        y_values_1.push_back(f(x_values_1.back()));
    }

    // 为定理3.58（二次B样条）定义节点
    int N2 = 10;
    for (int i = 1; i <= N2; ++i) {
        x_values_2.push_back(i - 11/2.0); 
        y_values_2.push_back(f(x_values_2.back()));
    }

    // 创建样条对象
    ClampedCubicSpline_B cubicSpline(3, x_values_1, y_values_1, 10.0/(26*26), -10.0/(26*26)); 
    QuadraticSpline_B quadraticSpline(2, x_values_2, y_values_2, 1.0/26, 1.0/26);

    // 绘制结果
    std::vector<double> x_plot;
    std::vector<double> y_cubic, y_quadratic, y_exact;

    double h = 10.0 /100;
    for (int i = 0; i < 100; ++i) {
        x_plot.push_back(-5 + i * h);
        y_cubic.push_back(cubicSpline.evaluate(x_plot.back()));
        y_quadratic.push_back(quadraticSpline.evaluate(x_plot.back()));
        y_exact.push_back(f(x_plot.back()));
    }

    // 使用 matplotlib-cpp 绘图
    // plt::named_plot( "Cubic Spline", x_plot, y_cubic, "r-");
    // plt::named_plot("Quadratic Spline", x_plot, y_quadratic, "b--");
    // plt::named_plot("Exact Function", x_plot, y_exact, "g-");
    // plt::legend();
    // plt::xlabel("x");
    // plt::ylabel("f(x)");
    // plt::title("Cubic vs Quadratic B-Spline Interpolation");
    // plt::save("F:/Codefield/numerical_analysis/Programming/Chapter3/figure/PC.png");
    // plt::show();

    // 定义插值误差计算的点
    vector<double> points = {-3.5, -3, -0.5, 0, 0.5, 3, 3.5};
    vector<double> errors_cubic, errors_quadratic;

    // 计算每个点的插值误差
    cout << "Interpolation Errors:" << endl;
    cout << "x\t\tkCubic\t\tkQuadratic" << endl;
    for (double point : points) {
        double error_cubic = calculateError(point, cubicSpline);
        double error_quadratic = calculateError(point, quadraticSpline);
        
        cout << point << "\t\t" << error_cubic << "\t\t" << error_quadratic << endl;
    }
    return 0;
}